home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / mshplanr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  9.3 KB  |  312 lines
  1. /******************************************************************************
  2. * MshPlanr.c - Test colinearity of control meshes/polygons.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Sep. 91.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. /*****************************************************************************
  10. * Computes a distance between two control points at given indices Index?.    *
  11. * Only E2 or E3 point types are supported.                     *
  12. *****************************************************************************/
  13. CagdRType CagdDistanceTwoCtlPts(CagdRType **Points, int Index1, int Index2,
  14.                             CagdPointType PType)
  15. {
  16.     switch (PType) {
  17.     case CAGD_PT_E2_TYPE:
  18.         return sqrt(SQR(Points[1][Index1] - Points[1][Index2]) +
  19.             SQR(Points[2][Index1] - Points[2][Index2]));
  20.     case CAGD_PT_E3_TYPE:
  21.         return sqrt(SQR(Points[1][Index1] - Points[1][Index2]) +
  22.             SQR(Points[2][Index1] - Points[2][Index2]) +
  23.             SQR(Points[3][Index1] - Points[3][Index2]));
  24.     default:
  25.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  26.         return 0.0;
  27.     }
  28. }
  29.  
  30. /*****************************************************************************
  31. * Fits a plane into the four points from Points indices Index?. Points may   *
  32. * be either E2 or E3 only.                             *
  33. *   Returns 0.0 if failed to fit a plane, otherwise a measure on the size of *
  34. * the mesh data (distance between points) isreturned.                 *
  35. *****************************************************************************/
  36. CagdRType CagdFitPlaneThruCtlPts(CagdPlaneStruct *Plane, CagdPointType PType,
  37.                      CagdRType **Points,
  38.                      int Index1, int Index2, int Index3, int Index4)
  39. {
  40.     int i, j, Indices[4];
  41.     CagdRType SizeMeasure,
  42.     MaxDist = 0.0;
  43.     CagdVecStruct V1, V2, V3;
  44.  
  45.     Indices[0] = Index1;
  46.     Indices[1] = Index2;
  47.     Indices[2] = Index3;
  48.     Indices[3] = Index4;
  49.  
  50.     /* Find out the pair of vertices most seperated: */
  51.     for (i = 0; i < 4; i++)
  52.     for (j = i + 1; j < 4; j++) {
  53.         CagdRType
  54.         Dist = CagdDistanceTwoCtlPts(Points, Indices[i], Indices[j],
  55.                                      PType);
  56.         if (Dist > MaxDist) {
  57.         MaxDist = Dist;
  58.         Index1 = i;
  59.         Index2 = j;
  60.         }
  61.     }
  62.  
  63.     if (MaxDist < EPSILON)
  64.     return 0.0;
  65.     SizeMeasure = MaxDist;
  66.     MaxDist = 0.0;
  67.  
  68.     /* Find a third most seperated than the selected two. */
  69.     for (i = 0; i < 4; i++)
  70.     if (i != Index1 && i != Index2) {
  71.         CagdRType
  72.         Dist1 = CagdDistanceTwoCtlPts(Points, Indices[Index1],
  73.                           Indices[j], PType),
  74.         Dist2 = CagdDistanceTwoCtlPts(Points, Indices[Index2],
  75.                           Indices[j], PType),
  76.         Dist = MIN(Dist1, Dist2);
  77.  
  78.         if (Dist > MaxDist) {
  79.         MaxDist = Dist;
  80.         Index3 = i;
  81.         }
  82.     }
  83.     if (MaxDist < EPSILON)
  84.     return 0.0;
  85.  
  86.     /* Now we have the 3 most seperated vertices to fit a plane thru. */
  87.  
  88.     switch (PType) {
  89.     case CAGD_PT_E2_TYPE:
  90.         /* This is trivial since the plane must be Z = 0. */
  91.         Plane -> Plane[0] = 0.0;
  92.         Plane -> Plane[1] = 0.0;
  93.         Plane -> Plane[2] = 1.0;
  94.         Plane -> Plane[3] = 0.0;
  95.         break;
  96.     case CAGD_PT_E3_TYPE:
  97.         /* Compute normal as a cross product of two vectors thru pts. */
  98.         for (i = 0; i < 3; i++) {
  99.         j = i + 1;
  100.         V1.Vec[i] = Points[j][Index2] - Points[j][Index1];
  101.         V2.Vec[i] = Points[j][Index3] - Points[j][Index2];
  102.         }
  103.         CROSS_PROD(V3.Vec, V1.Vec, V2.Vec);
  104.         CAGD_NORMALIZE_VECTOR(V3);
  105.         for (i = 0; i < 3; i++)
  106.         Plane -> Plane[i] = V3.Vec[i];
  107.  
  108.         Plane -> Plane[3] = (-(V3.Vec[0] * Points[1][Index1] +
  109.                    V3.Vec[1] * Points[2][Index1] +
  110.                    V3.Vec[2] * Points[3][Index1]));
  111.         break;
  112.     default:
  113.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  114.         break;
  115.     }
  116.  
  117.     return SizeMeasure;
  118. }
  119.  
  120. /*****************************************************************************
  121. * Computes and returns distance between point Index and given plane which is *
  122. * assumed to be normalized, so that the A B C D plane normal has unit len..  *
  123. * Also assumed the Points are non rational with MaxDim dimension.         *
  124. *****************************************************************************/
  125. CagdRType CagdDistPtPlane(CagdPlaneStruct *Plane, CagdRType **Points,
  126.                               int Index, int MaxDim)
  127. {
  128.     int i;
  129.     CagdRType
  130.     R = Plane -> Plane[3];
  131.  
  132.     for (i = 0; i < MaxDim; i++)
  133.     R += Plane -> Plane[i] * Points[i + 1][Index];
  134.  
  135.     return fabs(R);
  136. }
  137.  
  138. /*****************************************************************************
  139. * Computes and returns distance between point Index and the line from first  *
  140. * point in direction LineDir (usually the line direction to last the point). *
  141. * LineDir is assumed to be unit length normalized vector.             *
  142. *****************************************************************************/
  143. static CagdRType CagdDistPtLine(CagdVecStruct *LineDir, CagdRType **Points,
  144.                         int Index, int MaxDim)
  145. {
  146.     int i;
  147.     CagdRType R;
  148.     CagdVecStruct V1, V2;
  149.  
  150.     for (i = 0; i < MaxDim; i++)
  151.     V1.Vec[i] = Points[i+1][Index] - Points[i+1][0];
  152.     if (MaxDim < 3)
  153.     V1.Vec[2] = 0.0;
  154.  
  155.     /* Let V1 be the vector from point zero to point index. Then the distance */
  156.     /* from point Index the the line is: | (LineDir . V1) LineDir - V1 |.     */
  157.     V2 = *LineDir;
  158.     R = DOT_PROD(V1.Vec, V2.Vec);
  159.     CAGD_MULT_VECTOR(V2, R);
  160.     CAGD_SUB_VECTOR(V2, V1);
  161.  
  162.     return CAGD_LEN_VECTOR(V2);
  163. }
  164.  
  165. /*****************************************************************************
  166. * Test polygon colinearity by testing the distance of interior control       *
  167. * points from the line connecting the two polygon end points.             *
  168. *   Returns a relative ratio of deviation from line relative to its length.  *
  169. *   Zero means all points are colinear.                         *
  170. *   If two end points are same ( no line can be fit ) INFINITY is returned.  *
  171. *****************************************************************************/
  172. CagdRType CagdEstimateCrvColinearity(CagdCrvStruct *Crv)
  173. {
  174.     int i,
  175.     MaxDim = 3,
  176.     Length = Crv -> Length,
  177.     Length1 = Length - 1;
  178.     CagdRType LineLen,
  179.     MaxDist = 0.0,
  180.     **Points = Crv -> Points;
  181.     CagdCrvStruct
  182.     *CoercedCrv = NULL;
  183.     CagdPointType
  184.     PType = Crv ->PType;
  185.     CagdVecStruct LineDir;
  186.  
  187.     switch (PType) {          /* Convert the rational cases to non rational. */
  188.     case CAGD_PT_P2_TYPE:
  189.         CoercedCrv = CagdCoerceCrvTo(Crv, CAGD_PT_E2_TYPE);
  190.         Points = CoercedCrv -> Points;
  191.         PType = CoercedCrv -> PType;
  192.         break;
  193.     case CAGD_PT_P3_TYPE:
  194.         CoercedCrv = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE);
  195.         Points = CoercedCrv -> Points;
  196.         PType = CoercedCrv -> PType;
  197.         break;
  198.     default:
  199.         break;
  200.     }
  201.  
  202.     switch (PType) {
  203.     case CAGD_PT_E2_TYPE:
  204.         LineDir.Vec[0] = Points[1][Length1] - Points[1][0];
  205.         LineDir.Vec[1] = Points[2][Length1] - Points[2][0];
  206.         LineDir.Vec[2] = 0.0;
  207.         MaxDim = 2;
  208.         break;
  209.     case CAGD_PT_E3_TYPE:
  210.         LineDir.Vec[0] = Points[1][Length1] - Points[1][0];
  211.         LineDir.Vec[1] = Points[2][Length1] - Points[2][0];
  212.         LineDir.Vec[2] = Points[3][Length1] - Points[3][0];
  213.         break;
  214.     default:
  215.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  216.         break;
  217.     }
  218.  
  219.     LineLen = CAGD_LEN_VECTOR(LineDir);
  220.     if (LineLen < EPSILON) {
  221.     if (CoercedCrv != NULL)
  222.         CagdCrvFree(CoercedCrv);
  223.     return INFINITY;
  224.     }
  225.  
  226.     CAGD_DIV_VECTOR(LineDir, LineLen);
  227.  
  228.     for (i = 1; i < Length1; i++) {
  229.     CagdRType
  230.         Dist = CagdDistPtLine(&LineDir, Points, i, MaxDim);
  231.  
  232.     if (Dist > MaxDist)
  233.         MaxDist = Dist;
  234.     }
  235.  
  236.     if (CoercedCrv != NULL)
  237.     CagdCrvFree(CoercedCrv);
  238.  
  239.     return MaxDist / LineLen;
  240. }
  241.  
  242. /*****************************************************************************
  243. * Test mesh colinearity by testing the distance of interior points from the  *
  244. * plane thru 3 corner points.                             *
  245. *   Returns a relative ratio of deviation from plane relative to its size.   *
  246. *   Zero means all points are coplanar.                         *
  247. *   If end points are same ( no plane can be fit ) INFINITY is returned.     *
  248. *****************************************************************************/
  249. CagdRType CagdEstimateSrfPlanarity(CagdSrfStruct *Srf)
  250. {
  251.     int i,
  252.     ULength = Srf -> ULength,
  253.     ULength1 = ULength - 1,
  254.     VLength = Srf -> VLength,
  255.     VLength1 = VLength - 1;
  256.     CagdRType
  257.     PlaneSize = 0.0,
  258.     MaxDist = 0.0,
  259.     **Points = Srf -> Points;
  260.     CagdSrfStruct
  261.     *CoercedSrf = NULL;
  262.     CagdPointType
  263.     PType = Srf ->PType;
  264.     CagdPlaneStruct Plane;
  265.  
  266.     switch (PType) {          /* Convert the rational cases to non rational. */
  267.     case CAGD_PT_P2_TYPE:
  268.     case CAGD_PT_E2_TYPE:
  269.         /* Trivial case - it is planar surface. */
  270.         return 0.0;
  271.     case CAGD_PT_P3_TYPE:
  272.         CoercedSrf = CagdCoerceSrfTo(Srf, CAGD_PT_E3_TYPE);
  273.         Points = CoercedSrf -> Points;
  274.         PType = CoercedSrf -> PType;
  275.         break;
  276.     default:
  277.         break;
  278.     }
  279.  
  280.     switch (PType) {
  281.     case CAGD_PT_E3_TYPE:
  282.         PlaneSize = CagdFitPlaneThruCtlPts(&Plane, PType, Points,
  283.                     CAGD_MESH_UV(Srf, 0,        0),
  284.                     CAGD_MESH_UV(Srf, 0,        VLength1),
  285.                     CAGD_MESH_UV(Srf, ULength1, 0),
  286.                     CAGD_MESH_UV(Srf, ULength1, VLength1));
  287.         break;
  288.     default:
  289.         FATAL_ERROR(CAGD_ERR_UNSUPPORT_PT);
  290.         break;
  291.     }
  292.  
  293.     if (PlaneSize < EPSILON) {
  294.     if (CoercedSrf != NULL)
  295.         CagdSrfFree(CoercedSrf);
  296.     return INFINITY;
  297.     }
  298.  
  299.     for (i = ULength *VLength; i > 0; i--) {
  300.     CagdRType
  301.         Dist = CagdDistPtPlane(&Plane, Points, i, 3);
  302.  
  303.     if (Dist > MaxDist)
  304.         MaxDist = Dist;
  305.     }
  306.  
  307.     if (CoercedSrf != NULL)
  308.     CagdSrfFree(CoercedSrf);
  309.  
  310.     return MaxDist / PlaneSize;
  311. }
  312.